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Abstract 

This article is devoted to the study of an incompressible viscous flow of a fluid 
partly enclosed in a cylindrical container with an open top surface and driven by 
the constant rotation of the bottom wall. Such type of flows belongs to a group of 
recirculating lid-driven cavity flows with geometrical axisymmetry and of the pre- 
scribed boundary conditions of Dirichlet type — no-slip on the cavity walls. The top 
surface of the cylindrical cavity is left open with an imposed stress-free boundary 
condition, while a no-slip condition with a prescribed rotational velocity is imposed 
on the bottom wall. The Reynolds regime corresponds to transitional flows with 
some incursions in the fully laminar regime. The approach taken here revealed new 
flow states that were investigated based on a fully three-dimensional solution of 
the Navier-Stokes equations for the free-surface cylindrical swirling flow, without 
resorting to any symmetry property unlike all other results available in the litera- 
ture. Theses solutions are obtained through direct numerical simulations based on 
a Legendre spectral element method. 
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1 Introduction 



Besides the differences in terms of geometry, the hd-driven cubical cavity flow [1] and the cylin- 
drical swirling flow investigated in this article, present similar features typical of shear-driven 
recirculating flows such as intense wall-jets, shear layers in the vicinity of the driven wall, and 
secondary recirculating flows, all of which are very dependent on the flow parameters. Never- 
theless, the geometry — cubical on one hand and cylindrical on the other hand — dramatically 
influences the nature and structure of these secondary flows: corner eddies for the cubical cavity 
and recirculation bubbles or vortex breakdown in the cylindrical case. 

1.1 General considerations 

Following the pioneering work of Bogatyrev & Gorin [2] and Koseff & Street [3,4], it was shown 
that contrary to intuition, the lid-driven cubical cavity flow is essentially three-dimensional, 
even when considering large aspect ratio. It is only recently that the three-dimensionality of 
the lid-driven cylindrical cavity flow was confirmed numerically by Blackburn & Lopez [5, 6] 
after it was suggested but not fully proved experimentally by S0rensen [7], Spohn et al. [8], 
Sotiropoulos & Ventikos [9], and Pereira & Sousa [10]. In 2001, Sotiropoulos & Ventikos [11] 
gave full experimental evidence of the three-dimensional character of the flow with the onset of 
non-axisymmetric modes. The three-dimensional nature of these driven cavity flows therefore 
appears as a general characteristic of internal recirculating shear-driven flows. 

In the sequel, we will only consider the cylindrical lid-driven cavity flow also referred to as 
"swirling" flow without any additional precision. The flrst experiments by Vogcl [12] and later 
Ronnenberg [13] showed that Ekman suction and pumping, induced by the Ekman layers on 
the rotating and stationary disks, lead to the formation of a concentrated vortex core along the 
axis in the closed cavity case. The two dimensionless numbers characterizing this flow are the 
height-to-radius aspect ratio A = H/R and the Reynolds number Re = F^Q,/v, where H and 
R are the height and radius of the cylinder respectively, VL the constant angular velocity of the 
bottom end-wall, and v the kinematic viscosity of the Newtonian fluid. For speciflc values of 
the aspect ratio A, and above a critical Reynolds number, the vortex core breaks down in the 
form of one or more recirculation bubbles which are on-axis in the closed cavity case and on- 
or off-axis in the open cavity case. Owing to the enormous extent of work in the area of vortex 
breakdown (VB), (see reviews by Hall [14], Leibovich [15], Shtern & Hussain [16], KcrswcU [17], 
and Arndt [18]), we will only briefly recall the central features of VB. As defined by Leibovich 
in its review on the structure of vortex breakdown [15], the term "vortex breakdown" refers to 
a disturbance characterized by the formation of an internal stagnation point on the vortex axis, 
followed by reversed flow in a region of limited axial extent. Two forms of vortex breakdown 
predominate, one called "near-axisymmetric" (sometimes "axisymmetric" or "bubble- like" ) , and 
the other called "spiral" . 

The practical importance of vortex breakdown lies mainly in the fleld of aeronautics, where they 
can be observed over wings — mainly delta wings — with highly swept leading edges when the angle 
of incidence exceeds a critical value. Vortex breakdown can be a limiting factor on the operating 
altitude of slender- winged flying vehicles. Moreover, the occurrence of VB in the wake of a large 
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aircraft is relevant to the safety of flight in dense air-traffic, which is becoming more and more 
frequent with the constant increase in air-traffic over the years. VB is also important in other 
fields for example it has been observed in the swirling flows through nozzles and diffuscrs [19], 
and in the field of reactive flows, in combustion chambers. Besides the tremendous importance 
of VB in engineering applications, it is also a prototypical phenomenon allowing to elucidate the 
fundamental aspects of the bubble mode. 

1.2 The lid-driven cylindrical cavity flow 

The first comprehensive experimental study of the closed cylindrical container case was under- 
taken by Escudier [20], and Escudicr & Keller [21], who extended the earlier results of Vogel [12] 
and Ronnenberg [13] to obtain the first map of VB transitions with respect to the aspect ratio A 
and the Reynolds number. Escudier [20] revealed flow states with one, two or even three succes- 
sive breakdowns, as well as a transition from steadiness to unsteadiness. S0rensen [7] extended 
to a broader range of Reynolds number in the same experiment as Escudicr [20] for the closed 
container, and inferred that above a critical Reynolds number in the unsteady flow regime, the 
meridional flow becomes highly asymmetric. The flrst experimental study of the open cylindrical 
container case with a free surface on the top, was undertaken by Spohn et al. [22], who high- 
lighted the significant change in the structure, the occurrence and the location of the breakdown 
bubbles. In the steady closed cylinder case, Hourigan et al. [23] investigated the asymmetric 
spiraling effects along the cylinder axis prior to the flrst vortex breakdown. They argued that 
the observed asymmetry was purely an experimental artifact and not an evidence of the three- 
dimensional nature of the flow. Spohn et al. [8] were the first to investigate experimentally the 
origin of possible asymmetric features of the instabilities at their onset. The steady breakdown 
bubbles reported by Spohn et al. [8] showcase asymmetric features comparable to earlier mea- 
surements, and also to unsteady bubbles observed in circular diffusers by Faler & Leibovich [19]. 
As a matter of fact, the work of Spohn et al. [8] is really a pioneering work in the acceptance of 
the axisymmetry breaking, amongst fluid experimentalists, see Br0ns et al. [24]. It is noteworthy 
at this point, that the complex physics associated with these intricate phenomena occurring in 
closed/open rotating cylindrical container is still not clearly understood. 

Like for the lid-driven cubical cavity flow, and in relation with the simple geometry of the 
fiow, the rotating cylindrical cavity fiow has been extensively studied using direct numerical 
simulations. It is important to note that since the early seventies, the method of choice has 
consisted in solving the streamfunction-vorticity formulation of the axisymmetric incompressible 
Navier-Stokes equations. Without being exhaustive, the following list of references gives an 
overview of the numerical simulation of the closed lid-driven cylindrical flow over three decades: 
Pao [25], Lugt & Hausshng [26,27], Dijkstra & van Heijst [28], Lugt & Abboud [29], Neitzel 
[30], Daube & S0rensen [31], Lopez [32], Brown & Lopez [33], Lopez & Perry [34], S0rcnscn 
& Christensen [35], Watson & Neitzel [36], Gelfgat et al. [37,38], Tsitverblit & Kit [39], and 
Br0ns et al. [40]. All these works were able to reproduce with a reasonable accuracy, the basic 
features observed experimentally and reported earlier including the size, shape and number of 
recirculation bubbles. The onset of vortex breakdown was to some extent captured by several of 
these numerical simulations. Lopez [32], and Brown & Lopez [33] suggested a physical mechanism 
for the intricate phenomena observed. They prove the existence of a standing centrifugal wave, 
whose amplitude increases with the Reynolds number and which can create a stagnation point on 
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the cylinder axis, initiating the breakdown process. It is worth recaUing that the streamfunction- 
vorticity formulation is adequate and appropriate only for the study of flow dynamics preserving 

the property of axisymmctry. At the inception of any instability breaking the axisymmctry 
of the flow, a three-dimensional solution of the Navier-Stokes equations is required, thereby 
increasing considerably the complexity of the task. The last remark justifies the observed changes 
in terms of numerical modeling of Lopez' group and S0rensen's group, to allow them to investigate 
axisymmctry breaking in the closed cylinder case [5,6,41-43]. Therefore, three-dimensional flow 
structures have started being simulated more recently, see Gelfgat et al. [44], Sotiropoulos & 
Ventikos [11], Sotiropoulos et al. [45], Marques & Lopez [46], Blackburn & Lopez [5,6], Serre & 
Bontoux [47], Blackburn [48], and Lopez [49]. 

Apart from the canonical case with a single driving lid in rotation at a constant angular velocity, 
different variations of the problem have been extensively studied in the past years: e.g. cylinder 
with CO- and counter-rotating end-covers by Br0ns et al. [40], steady axisymmetric flow in an 
open cylindrical container with a partially rotating bottom wall by Piva & Meiburg [50] , vortex 
scenario and bubble generation in a cylindrical cavity with rotating top and bottom by Okulov 
et al. [51]. MuUin et al. [52] also included a rod at the axis to control the breakdown, and Pereira 
& Sousa [10] signiflcantly changed the conflguration by replacing the flat rotating bottom cover 
by a cone. As noted by Br0ns et al. [53], all these studies show a large set of flow structures 
which are quite sensitive to variations of external parameters. Mununga et al. [54] and Lo Jacono 
et al. [55] investigated different strategies for the control of vortex breakdown. 



1.3 Open swirling flow 

The focus in the present article is on the canonical problem of a cylinder with a rotating bottom 
end-wall but replacing the stationary sohd top end-wall by a free surface. The flow associated 
with this problem was flrst studied experimentally by Spohn et al. [8, 22] . They observed the 
influence of the top free surface — assuredly clean of surfactants — on the onset, structure, nature 
and number of recirculating bubbles. Their central observations arc that breakdown bubbles still 
appear, but are off-axis and may be attached to the free surface, depending on the aspect ratio 
A and the Reynolds number. Of course, such structures could not be observed in the closed case 
because of the no-slip condition imposed on the top wall. All the past simulations of free-surface 
swirling flows rely on the central assumptions that the free surface is flat and clean, which means 
that the Froude number is very small and that surface tension effects are negligible. With these 
assumptions, the flow is identical to the flow in the lower half part of a cylinder with two solid 
covers in co-rotation, i.e. rotating at the same angular velocity. Br0ns et al. reported a wide 
range of topologies of vortex breakdown bubbles in a bottom-driven cylinder with a free surface. 
Valentine & Jahnke [56], observed in their simulations the existence of one or two toroidal-like 
types of recirculation bubble having their stagnation lines attached to the free surface, depending 
on the value of the Reynolds number. Their study was complemented by the work of Lopez [57] 
for oscillating unsteady flows. Information relevant to the present problem with a free surface 
all indicate consistent flow behavior at small aspect ratio i.e., 0.5 < A < 1.0 in that stagnation 
occurs off-axis and associated secondary flow creates a toroidal recirculation bubble. Steady free- 
surface flows have been computed by Iwatsu [58,59] providing flow state classiflcations with new 
flow patterns not revealed in the previous studies. 
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1.4 Motivations and objectives 



The present study is motivated by several factors. Firstly, compared to the closed cylinder case, 
only some limited aspects of the open swirling flow have been investigated so far. The study of 
this intricate problem is relatively new and consequently the body of knowledge in some (A, Re)- 
parameter regions appears fairly limited. Secondly, most of the past studies involving numerical 
simulations of this free-surface swirling flow, used axisymmetric streamfunction-vorticity formu- 
lations: Br0ns et al. [53], Iwatsu [58,59], and Piva & Meiburg [50]. To our knowledge, the only 
fully three-dimensional numerical simulations is due to Lopez et al. [60] , who investigated mainly 
symmetry breaking issues. 

In the present article, new flow states are investigated based on a fully three-dimensional solution 
of the Navier-Stokes equations without the need to resort to symmetry properties by doubling the 
computational domain and enforcing co-rotation of both end-walls. To our knowledge, the present 
study provides the most general available results for this fiat-free-surface problem. Both, steady 
and unsteady flows are considered for different sets of governing parameters (A, Re). A Legendre 
spectral element method is used to provide an accurate solution of the governing equations, 
while the stress-free boundary condition is naturally enforced into the weak formulation of the 
problem. 

The mathematical model and the problem formulation are detailed in Sec. [21 while the original 
computational approach of this study is presented in Sec. [31 Subsequently, Sec. [H contains all 
the numerical results corresponding to various physical situations and flow states. Finally, the 
article ends with Sec. [5l providing summary and conclusions on the present work. 



2 Mathematical model and problem formulation 



2.1 Mathematical description of the problem 



The fluid enclosed in the cylindrical cavity is assumed to be incompressible, Newtonian with 
uniform density and temperature. The flow is governed by the Navier-Stokes equations 

dt dxj dxj 



duj_ 
dxj 



0, (2) 



where (7*j = —p5ij + 2vDij is the reduced Cauchy stress tensor of the fluid, p the static or reduced 
pressure, Dij the rate-of-deformation tensor, v the assuredly constant and uniform kinematic 
viscosity, and gi the components of the acceleration of gravity {gi = g2 = and g^ = —g). Inside 
the fluid domain denoted by V, no-slip boundary conditions are imposed on all cavity walls: 
the tubular side-wall, the bottom end-wall in steady rotation, see Fig. [H The mathematical 
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Fig. 1. Schematic of the geometry studied with the set of coordinates employed, 
expression of the no-slip condition on the tubular side-wall simply reads 

u{r = R, e, z, t) = v{r = R, 6, z, t) = w{r = R,e, z,t) = 0, < < 27r, < 2 < H. (3) 

The flow is driven by imposing a prescribed angular velocity distribution of the bottom end-wall, 
which transfers its kinetic energy to the fluid above. The details regarding the imposition of this 
Dirichlet boundary condition for the velocity field at the lid are discussed in the next section [2^21 
The top surface is left open and is modeled as a flat, fixed and clean free surface. The details 
regarding the imposition of this stress-free condition on the free surface are discussed in Sec. 12.31 
As mentioned in the Introduction, Sec. 11.11 two parameters that determine completely the flow 
state are the height-to-radius aspect ratio A = H/R and the Reynolds number Re = R'^Qq/u, 
based on the maximal prescribed angular velocity Qq of the bottom end-wall. 

In the sequel, the length, time, velocity, vorticity, helicity, streamfunction, (reduced) pressure 
and kinetic energy, and enstrophy are non-dimensionalized with respect to the reference scales 

Rj f2g , -Rfio) ^0; 



R^Ql, Ql respectively. 



2.2 Angular velocity distribution 



As already mentioned in Bouffanais et al. [61] for the study of the lid-driven cubical cavity flow, 
imposing a given angular velocity distribution on the bottom end-wall of a cavity is not an 
easy task numerically. Indeed, imposing a constant angular velocity profile leads to a singularity 
(discontinuous behavior in the velocity boundary conditions) at the circular edge between the 
bottom end-wall and the tubular side-wall, see Fig. [H Without adequate treatment, this discon- 
tinuous behavior will undermine the convergence and the accuracy of any numerical method in 
the vicinity of the lid. The same remedy as in the lid-driven cubical cavity problem in [61] is 
used here for the same reasons and with analogous justifications. A regularized angular velocity 
profile is employed by prescribing the following high-order polynomial expansion which vanishes 
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along its first derivatives on the circular edge 



n{r,9,z = 0,t) = Qo 



16 



(4) 



which leads to the following expressions in Cartesian coordinates of the components of the 
prescribed velocity field on the bottom end-wall 



u{x, y,z = 0,t) = u^{x, y,z = 0,t) = -yQo 

v{x, y, z = 0,t) = Uy{x, y, z = 0,t) = +xQo 
w{x, y,z = 0,t) = u^{x, y,z = 0,t) = 0, 



1 - + y^R 



l-[Jx'^ + y^R 



16 



16' 



(5) 

(6) 
(7) 



where x = {r,6 = 0) and y = {r,6 = 7r/2). This profile flattens very quickly near the circular 
edge {r/R = l,9,z/H = 0) while away from it, it grows rapidly to a constant value Qq of 
the angular velocity over a short distance. The highest polynomial degree of this distribution 
is 32. Such high-order polynomial expansions lead to steep velocity gradients in the vicinity 
of the circular edge of the bottom end-wall. The grid refinement, in terms of spectral element 
distribution near the disk will be presented in greater details in Sec. [31 One of the constraint in 
the grid design is to ensure the proper resolution of the lid velocity distribution by the spectral 
element decomposition. 



2.3 Free-surface modeling 



The analysis of a two-phase flow is based on the coupled hydrodynamics interactions between 
adjacent layers over a broad range of space and time scales. This analysis can be significantly 
simplified if the dynamics of the interface is almost entirely dependent — from the hydrodynamics 
and physico-chemistry viewpoints — on one phase, e.g. a liquid phase, and almost independent of 
the dynamics of the second phase, e.g. a gas phase. Based on this hypothesis, the surface is said 
to be free. Consequently, the two fluid phases can only exert constant normal stresses. Sarpkaya 
in his review entitled "Vorticity, free surface and surfactants" [62] gives a clear characterization 
and definition of a free surface: "Although, the exterior of a free surface is free from externally 
imposed shear, the interior is not necessarily free from the shear generated internally. In fact, 
surface deformations and contaminants give rise to surface-gradients and tangential stresses 
in the internal side of the bounding interface. From a mathematical viewpoint, a free surface 
means that the density and the viscosity of the upper fluid are zero and that the existence of a 
continuum above the interface is inconsequential. From a practical point of view, the free surface 
means that the dynamics of the continuum above the interface has negligible influence on the 
lower phase, i.e. a free surface is a simplifying approximation for an almost free surface." 

In the present study, the modeling of the interface between the fluids in the cylindrical cavity 
as a free surface is supplemented by an additional simplifying approximation: the free surface is 
supposed to remain flat and fixed all along the dynamic range of investigation. In general, the 
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dynamics of the free surface depends on the non-dimensional Froude number defined here as 



= (8) 

which measures the relative importance of the inertial effects compared to the stabilizing gravi- 
tational effects. Therefore, assuming a fiat free surface corresponds mathematically to a Froude 
number identically zero. As a consequence and in consistency with the latter assumption, the 
axial component of the velocity w = needs to vanish at the free surface z/H = 1 

w{x,y,z = H,t) =0, x^ + y^<R^, (9) 

thereby expressing the kinematic boundary condition at the free surface. This latter condition 
on the axial velocity w is to be supplemented with the stress-free condition at the free surface 

CTijfij = —p6ijhj + 2h'Dijnj = 0, (10) 

where n is the local outward unit vector normal at the free surface, which in the present particular 
situation is the unit normal vector in the z-direction. Consequently, the stress-free condition 
ffTU]) reduces to a*^ = a*^ = 0, i = 1,2,3 at z = H, and is explicitly stated as 

al, = 2uDr, = 0, (11) 

^2*3 = 2Z/D23 = 0, (12) 

a3*3 = -p + 2z/D33 = 0, (13) 
which under the zero-deformation condition simplifies to 

0, (14) 

0, (15) 

0. (16) 



du 


dw\ 
dx J 


du 


dz 


dz 


' dv 


dw\ 
dy J 


dv 


dz 


dz 




-p + 2z/ 


dw 
dz 



3 Computational approach 



3.1 Space discretization 



The Navier-Stokes equations ([I])-(l2]), supplemented with the boundary conditions ([3]), (l5])-(I7j), 
&, (fT^ - (fT6l) . constitute the set of governing equations for this free-surface swirling flow to be 
discretized and ultimately solved. The numerical method treats Eqs. ([l])-([2]) within the weak 
Galerkin formulation framework. The spatial discretization uses Lagrange-Legendre polynomial 
interpolants. The reader is referred to the monograph by Deville et al. [63] for full details. The 
velocity and pressure are expressed in the Pat — IPAr_2 functional spaces where Fjy is the set of 
polynomials of degree lower than N in each space direction. This spectral element method avoids 
the presence of spurious pressure modes as it was proved by Maday & Patera [64] and Maday et 
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al. [65]. The quadrature rules are based on a Gauss-Lobatto-Legendre (GLL) grid for the velocity 
nodes and a Gauss-Legendre grid (GL) for the pressure nodes. The spectral element grid used 
for all simulations is presented on Fig. O in the particular case A = 1. This mesh comprises 440 
spectral elements distributed into 10 cylindrical layers of different heights, but all made of the 
same distribution of 44 spectral elements, see Fig. [2] (right). In order to resolve the boundary 
layer along the tubular side-wall, the Ekman shear layer above the rotating bottom end-wall and 
the surface shear layer below the free surface, the spectral elements are unevenly distributed as 
can be seen in Figs. [2l The choice of polynomial degree in the three space directions, defining the 
inner GLL and GL grid into each spectral element is deferred to Sec. 13.31 The essential Dirichlet 
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Fig. 2. Typical meshes view used throughout this study. Left: half-meridional grid. Center: spectral-ele- 
ment grid in any plane normal to the z-direction. Right: Three-dimensional grid comprising 10 cylin- 
drical layers of nonuniform heights made of 44 spectral elements each. Case A = 1. 

boundary conditions — homogeneous for u on the tubular side- wall as expressed by Eq. ([3]), 
homogeneous for w = on the free surface as expressed by Eq. iQ, and non- homogeneous for 
u on the rotating bottom end-wall as expressed by Eqs. ©-([T]) — are embodied into the choice 
of test and trial functions chosen for the velocity field. 

The stress- free condition ( ITOi) at the free surface, further expressed by Eqs. (fT^-( fT6|l . appears 
as a homogeneous natural boundary condition in the weak Galerkin framework. This central 
point is discussed in detail and in a more general framework, in Bodard et al. [66] , where a non- 
homogenous natural boundary condition is accounted for in the weak formulation of the problem. 
Based on this previous analysis, the treatment of the present stress-free condition at the free 
surface appears straightforward and is automatically incorporated into the weak formulation of 
the problem. 

Borrowing the notation from Deville et al. [63], the semi-discrete Navier-Stokes equations re- 
sulting from space discretization are 

+ Cu + z/Au - D^p = 0, (17) 
at - 

-Du = 0. (18) 

The diagonal mass matrix M is composed of three blocks, namely the mass matrices M. The 
global vector u contains all the nodal velocity components while p is made of all nodal pressures. 
The matrices A, D^, D are the discrete Laplacian, gradient and divergence operators, respec- 
tively. The matrix operator C represents the action of the non-linear term written in convective 
form u ■ V, on the velocity field and depends on u itself. The semi-discrete equations constitute 
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a set of non-linear ordinary differential equations (fTTI) subject to the incompressibility condition 



3.2 Time integration 

The time discretization of the semi-discrete set of governing equations (fT71) - (fT8l) is the same as the 
one previously used in Bouffanais et al. [61] and Habisreutinger et al. [67]. We only briefly recall 
here the fundamentals of the method. The set of semi-discrete equations (fT71) - (|T8l) is discretized 
in time using finite-difference schemes in a decoupled approach. The computation of the linear 
Helmholtz problem — corresponding to the stiffness matrix A — is integrated based on an implicit 
backward differentiation formula of order 2, the nonlinear convective term — corresponding to 
the operator C — is integrated based on a relatively simple extrapolation method of order 2, 
introduced by Karniadakis et al. [68], see Bouffanais [69] for full details. 



3.3 Convergence tests 

In order to demonstrate the spatial and temporal convergence of the simulation method, time 
series data have been analyzed, while varying separately the time-step At and the polynomial 
degree of the GLL basis in each space direction, at the upper bound in Reynolds number 
Re = 6 000, and with A = H/R = 1. As no experimental, nor numerical reference results are 
available for the present problem, three integral and one local quantities have been computed 
and compared. These three integral quantities are the total kinetic energy Q, enstrophy E and 
helicity Ti. of the flow, which definitions are recalled 



The local quantity monitored is the axial velocity component w = at the point Ho of co- 
ordinates [r/R = 0,z/H = 0.75), located along the cylinder axis. The location choice of this 
monitoring point is motivated by the study of Piva & Meiburg [50] for a very similar configu- 
ration but at smaller Reynolds numbers. They show that in the vicinity of this point, the axial 
velocity component reaches a local maximum. Given the relatively high Reynolds number of our 
benchmark simulation, a quite long transient — approximately 500 time units in fig ^ units — is 
observed. Performing convergence tests involving a simulation time of the order of this tran- 
sient time would simply be prohibitive. Consequently, it was chosen to assess the convergence 
after only 50 time units of simulations, which corresponds to the appearance of the first vortex- 
breakdown recirculation bubble in the fluid domain. First, the spatial convergence is assessed by 
varying the polynomial degree in the range 6 < < 10, while keeping the time-step values con- 
stant At = 0.0025. Results are reported in Tabled! and suggest that the spatial convergence is 
achieved using a polynomial degree A^ = 8 in all three space directions. This value is used for all 
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(20) 
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the direct numerical simulations presented in the sequel, except for one single case corresponding 
to the steady laminar case (A = 1, Re = 900), for which iV = 7 is chosen. 
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Table 1 

Spatial convergence analysis for the case (Re = 6 000, A = 1) with At = 0.0025 J^^^ Q in R^{mo)^ 
units, E in R^^l^ units, 7i in R^^Iq units, and w in Rfto units. Instant t = 50Qq^. 

Finally, the temporal convergence is assessed by varying the time-step, while keeping the poly- 
nomial degree constant at the value = 8, in agreement with the earlier spatial convergence 
analysis. Results are reported in Table [2] below, and suggest that the temporal convergence is 
achieved using a time-step At = 0.0025. At a reduced Reynolds number compared to the one 
employed for this convergence analysis Re = 6 000, greater values of the time-step have been 
chosen in relation with the more laminar nature of the flow without affecting the convergence of 
the simulations. 

At Q E n wiUo) 

0.0050 2.08574e-02 2.14302e+01 1.36342e-01 3.49221e-03 

0.0035 2.11896e-02 2.25311e+01 1.48303e-01 3.61923e-03 

0.0025 2.19036e-02 2.37953e+01 1.66448e-01 3.81373e-03 

0.0010 2.19034e-02 2.37960e+01 1.66446e-01 3.81379e-03 

Table 2 

Temporal convergence analysis for the case (Re = 6000, A = 1) with N = 8. At in 0,^^ units, Q in 
i?^(i?r2o)^ units, E in R^Qq units, 7i in i?^r2g units, and w in R^Iq units. Instant t = 50^1^^. 



4 Numerical simulations and results 

4-.1 General physical characteristics of the flow 

The central characteristics of the flow in a closed cylindrical container with a bottom rotating 
end-wall is a large recirculation of the fluid. The features of the intense shear layer induced 
by the rotation of the bottom wall can be obtained from the analogy with the analysis by 
von Karman for the flow generated by a spinning plate in an unbounded fluid domain, see the 
review by Zandbergen & Dijkstra [70] for full details. The rotation of the bottom wall has a 
suction effect on the fluid in the near- axis region and a pumping effect, while accelerating the 
fluid radially outwards in an Ekman shear layer of thickness 0(Re~^'^^). In the framework of 
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our problem, this Ekman layer is bounded by the tubular cylinder side-wall, which forces the 
recirculation of the fluid in the upward direction along the side- wall, and towards the top wall. 
As the fluid approaches the curved corner, the radial velocity contribution to the kinetic energy 
is progressively transformed into an axial velocity contribution to the kinetic energy. It has been 
observed that the fluid turns and subsequently spirals upward along the tubular side-wall. 

As mentioned in Sec. 11.31 replacing the fixed top solid wall with a free surface changes significantly 
the physics of the flow and the recirculation mechanisms. In absence of tangential stresses at the 
free surface, the boundary layer is replaced by a surface layer in the sense of Shen et al. [71,72]. 
In addition, the inward spiraling fluid elements conserve their angular momentum at the free 
surface. When the related centrifugal force is large enough to balance the radial pressure gradient, 
the flow separates from the free surface and leads to the generation of a vortex breakdown bubble. 
The most striking difference between the flow patterns observed in the present open cylinder case 
and compared to the close cylinder one are the possible appearances of recirculation bubbles, 
which are generally attached to the free surface. Such flow patterns are simply impossible in 
presence of the no-slip condition imposed on the top wall in the closed cylinder case. More 
precisely, Iwatsu [58] determined 24 different flow states in the steady regime according to the 
meridional streamline patterns observed. Spohn et al. [8] summarized those flow states in a 
simplified bifurcation diagram in the (A, Re) space. Based on his extensive and comprehensive 
study, Iwatsu [58] came out with a more detailed and complex bifurcation diagram. 

In subsequent studies, the stability of those steady axisymmetric flows were investigated. Young 
et al. [73], Miraghaie et al. [74], Lopez et al. [60], and Lopez & Marques [75] observed unstable 
azimuthal modes which are triggered at different values of the Reynolds number depending on the 
nature "shallow" (A < 1) or "deep" (A > 1) of the system. Valentine & Jahnke [56], Lopez [57], 
and Br0ns et al. [53] associated the axisymmetry breakage to instability modes. These modes 
appear following a Hopf bifurcation which generally occurs at relatively high Reynolds number. 

4-2 Cases studied 

As mentioned in the previous sections, the physics of these free-surface swirling flows depends 
critically on the Reynolds number. Nevertheless, the height-to-radius aspect ratio A also has 
considerable impact on the observed nature of the flow. Very often, situations corresponding to 
extreme values of A have been studied, as they generally lead to simplified flow mechanisms. For 
instance, shallow systems (A < 1) are often referred to as "rotor-stator" configurations, in which 
the fluid is almost in a state of solid-body rotation. On the other hand, deep systems associated 
with large values of A, generate recirculation bubbles away from the free surface and generally 
on the cylinder axis. Consequently, systems corresponding to values of A close to the unity are 
intermediate in the sense that the physics of the flow observed is a complex combination of the 
general features reported for the shallow and deep systems. 

The details related to the five cases considered in this article are summarized in Table [31 The 
primary focus is on the flow defined by Re = 6 000 and A = 1, and corresponding to case (c). 
The value of the Reynolds number is intentionally set to a high value compared to previous 
studies — the highest to our knowledge, in order to obtain fields of a relative intensity at the free 
surface. The choice of the value of A follows the earlier comment on systems being intermediate 
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Case 


Re 


A = H/R 


Time-step 


Time evolution 


Vortex breakdown 


(a) 


900 


1.0 


0.0050 


steady 


one attached bubble 


(b) 


1500 


1.0 


0.0050 


steady 


one attached toroidal bubble 


(c) 


6 000 


1.0 


0.0025 


unsteady 


comolex dynamics 


id) 


2 000 


1/3 


0.0040 


steady 


two long attached bubbles 




2 000 


3.0 


0.0040 


steady 


one short detached bubble 



Table 3 

Parameters and characteristics of the cases considered. The time-step At is expressed in JIq ^ units. 
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Fig. 3. Time history of the volume integral of the kinetic energy Q of the flow, in R^Q,q units for cases 
(a)-(e). 

between shallow and deep. This central case (A = l,Re = 6 000) is supplemented with four 
secondary cases described in Table [31 The study of those secondary flows is of prime importance 
for the understanding of the complex dynamics of the primary case (A = 1, Re = 6 000). 

In terms of initial conditions, the steady rotation is impulsively started from a quiescent fluid 
state for all cases presented in the sequel. At this point, it is worth noting the timescales of the 
evolution of these flows. Figure [3] displays the time history of the volume integral of the kinetic 
energy of the flow Q. For all cases except case (c), the flow reaches a steady state after a given 
time scale, which is, as expected, shorter for shallow systems. Case (c) leads to an unsteady flow 
which does not display any oscillatory evolution. The value of the Reynolds number for this case 
is large enough to produce a non-trivial evolution of the recirculation zones as will be seen in 
the sequel. 

The time histories of the volume integral of the kinetic energy Q for the five cases (a) — (e) can 
be compared to the ones, reported in Bouffanais [69] for the closed swirling flow problem with 
A = 2.5, see Fig B.4. It should be noted that for a fixed value of A, the total kinetic energy Q 
of the flow decreases with the Reynolds number for the closed cylinder case, while it increases 
in the open cylinder case. This decreasing trend for Q in the closed cylinder case can easily be 
resolved by transposing the analysis given by Leriche & Gavrilakis [76] in their study of the 
closed lid-driven cubical cavity flow. Leriche & Gavrilakis argue that the most significant part 
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of the kinetic energy of the flow is contained in the viscous layer developing on the driving wall. 
Consequently, the total energy varies like the energy contained in this viscous layer which can 
approximately be expressed as f/g VRe~^^^, where Uq is the characteristic velocity of the driving 
wall and V the volume of the cavity. Such argument and estimate can easily be transposed for the 
closed swirling flows, and explains the decreasing trend for Q with respect to Re. Furthermore 
this argument is confirmed by the measurements of the kinetic energy Q{Li) of the cylindrical 
layer of fluid Li located right above the spinning disk and of height O.OlSif, reported in Table |H 
This thin layer of fluid which only represents 1.5% of the total volume of fluid, contributes for 
approximately 10% to the total kinetic energy of the flow. Concurrently, its contribution to the 
total kinetic energy of the flow decreases with the Reynolds number. 



Re 




Q{Li) 




E{Li) 






900 


7, 


.42243e-03 


3. 


.04128e+00 


3. 


,99817e-01 


1500 


7, 


,05011e-03 


4. 


.46885e+00 


4. 


,03147e-01 


6 000 


5, 


,71943e-03 


1, 


. 18305e+01 


4, 


,25622e-01 


Re 




Q{Lio) 








Ez{Lio) 


900 


3, 


.88037e-04 


6. 


. 14708e-03 


5. 


,99848e-03 


1500 


5, 


,47285e-04 


9. 


.04519e-03 


8. 


,25633e-03 


6 000 


1, 


.21966e-03 


6. 


.43344e-02 


2, 


,00737e-02 



Table 4 

Measurements of the kinetic energy Q, enstrophy E, and enstrophy associated solely with the axial 
vorticity component for the cylindrical layer Li comprised between z = and z = 0.015/7, and for 
the cylindrical layer Liq comprised between z = 0.98H and z = H. Q in R^{RQq)'^ units, E and E^ in 
R^nl units. Instant t = 6001^0 ^ 

On the contrary, a reverse trend is observed for the variations of Q with respect to Re in the 
open swirling flow. It therefore requires another physical justification. Nevertheless, the previous 
energetic argument associated with the viscous layer still holds for the viscous layer above the 
spinning disk and near the tubular side-wall in the open cylinder swirling flow. Below the surface 
at z = H, the viscous layer in the closed cylinder case is replaced by an intense shear layer. As 
mentioned in Sec. 12.31 in the present flat-free- surface problem, the axial vorticity is the only 
component of the vorticity field which is not vanishing at the free surface and it provides a 
measure of the internal shear at the free surface. The part of the enstrophy associated with 
the axial vorticity component is measured in the cylindrical layer of fluid Liq, of height 0.02H 
and located below the free surface. Results are reported in Table H] and clearly show a significant 
increase of Ez with respect to Re in the layer Lio, while it is almost constant in Li. These results 
allow us to infer that the shear layer below the free surface becomes more and more intense 
and energetic — see Q(Lio) — when increasing the Reynolds number. But this observed energetic 
trend of the free-surface shear layer is not the only factor responsible for counterbalancing the 
decreasing trend of the viscous layers. The internal structure of the free-surface layer is itself 
physicahy different as reported by Shen et al. [71,72]. The so-called surface layer corresponds to a 
thin region adjacent to the free surface characterized by fast variations of the tangential vorticity 
components. This surface layer is caused by the dynamic zero-stress boundary conditions ffTOl) 
at the free surface and lies inside a thicker blockage (or "source") layer, which is due to the 
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kinematic boundary condition (Q at the free surface. The importance of the outer blockage 
layer is manifested mainly in the redistribution of the kinetic energy, i.e. in the increase of 
the kinetic energy of the tangential velocity components at the expense of the kinetic energy 
of the axial velocity component. This point is to be further discussed in Sec. 14.3.21 where a 
comprehensive comparison of the flow below the surface z = H for the case (c) with a free 
surface and the equivalent closed case is given. 

4.3 Physical description of flow states 

4.3.1 Steady flows 

As a first step, we present the two steady flows for A = 1 at Re = 900 and 1 500 corresponding to 
cases (a) and (b) respectively. Figure H] displays the streamlines of these flows into any meridional 
plane. Both of these flows present a large axisymmetric vortex breakdown bubble attached to 
the free surface in agreement with the experimental results from Spohn et al. [8] summarized in 
their bifurcation diagram. These recirculation zones are characteristic of these swirling flows due 
to the conjugate effects of the centrifugal force and the overturning flow induced by the presence 
of the tubular side-wall. The central difference between the low-Reynolds-number cases (a) and 
(6) is the shape of the recirculation, which becomes toroidal after leaving the axis when the 
Reynolds number is increased from 900 up to 1 500. These results can be further validated by 
comparing them to the experimental results (dye visualizations) obtained by Piva for Re = 1 120. 
Finally, case (a) = (A = 1, Re = 900) has also been computed by Piva & Meiburg [50]. 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

Fig. 4. Contours of streamlines in a meridional plane, case A = 1. Left: case Re = 900; Right: case 
Re = 1 500. The 30 contours are non-uniformly spaced for visualization purposes, 20 equally-spaced 
negative contours and 10 equally-spaced positive contours for Re = 900 and Re = 1 500. 

As a second step, the contours of the radial, azimuthal and axial velocity components in any 
meridional plane are given in Fig. [HI These data are supplemented with the contours of the 
axial component of the angular momentum F = ru0 still in Fig. [HI extreme right column. 
The interest for F lies in the fact that it plays the role of a streamfunction for the part of the 
velocity field comprised in any meridional plane, see Bragg & Hawthorne [77] and Keller [78] 
for full details. Therefore, the contours of F deliver the intersection of vortex surfaces with the 
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corresponding meridional plane where they are drawn, and as such provide us with the local 
direction of the meridional vorticity field. One can notice from the velocity components and 
axial angular momentum component that the meridional structure of these fiows is far from 
being trivial. It consists of an intense boundary layer above the spinning bottom end-wall that is 
turned into the interior by the presence of the tubular side-wall, forming a shear layer having a 
jet-like velocity profile in the azimuthal direction. The contour lines of the axial component of the 
angular momentum shown in Fig. [T3] (extreme right column) simply represent the vortex lines, 
which all emanate from the rotating disk; the structure of the shear layer is apparent. It is worth 
noting here that the vortex lines distribution at their origin varies like r^. As a consequence of 
the regularized profile of angular velocity of the rotating disk — see Sec. \2.2\ this distribution in 
is slightly affected in the vicinity of r = R. This regularization of this profile has the advantage 
of preventing the appearance of vortex lines terminating at the circular corner {r = R,z = 0). 
The overturning nature of these flows is also apparent in the vicinity of the tubular side-wall, 
which is the vortex surface corresponding to F = 0, together with the cylinder axis. As non-zero 
azimuthal velocities are possible at the free surface, vortex lines emanating from the rotating 
end-wall have the option of terminating orthogonally to the free surface. This observation is 
one of the major difference with the closed cylinder swirling flow where all vortex lines have 
to terminate in the corner. Furthermore, the termination of vortex lines at the free surface is 
responsible for the possibility of having vortex breakdown bubbles being attached to the free 
surface as observed in Fig. HI A careful analysis of Fig. [H] (extreme right column) reveals that 
one vortex line marks the limit between an inner region comprising only vortex lines terminating 
at the free surface, and an outer region, where they terminate near the circular corner, like in 
the closed cylinder case. In summary, it appears that the main effect of this overturning flow is 
to bring high- angular-momentum fluid towards the cylinder axis. 

The results for case (a) = (A = l,Re = 900) presented in Fig. H] (left) and Fig. [H] (top row) 
show a qualitative good agreement with the numerical results of Piva & Meiburg [50]. The 
differences related to the features of the recirculation bubbles for cases (a) and (6) have been 
discussed earlier. One can notice in Fig. [H] (two flrst left columns), that the thicknesses of the 
intense radial velocity layer as well as the axial wall jet are reduced when the Reynolds number 
is increased from 900 to 1 500, as expected. The contours of the axial velocity component reveal 
that the downward-directed flow induced by the suction effect of the Ekman layer, is more 
intense at higher Reynolds number. In addition, the region of the flow where w = Uz has a 
negative extremum tends to move closer to the free surface when increasing Re. Regarding the 
vortex lines shown in the extreme right column, their bending towards the cylinder axis is more 
pronounced at the higher Reynolds number of 1 500. In relation with the previous analysis, this 
latter observation highlights the fact that more high-angular-momentum fluid is brought towards 
the axis when increasing Re. 

As mentioned in Sec. 14.21 our primary interest lies in case (c) = (A = l,Re = 6 000), thereby 
justifying the study of cases (a) and (6), having the same aspect ratio A but corresponding to 
laminar cases. Nevertheless, the study of cases {d) and (e), which both correspond to "extreme" 
cases in terms of height-to-radius aspect ratio, illustrate some essential features of the open 
swirling flow. In case (c) some of these features may prevail only in specific regions of the fiow as 
it corresponds to an intermediate case between a shallow system characterized by case {d) and 
a deep system characterized by case (e). These features cirG clS follows: 
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- solid-body rotation of the inner core region, predominantly for small A; 

- radial jet of angular momentum at the free surface; 

- lateral jet-like shear layer along the tubular side-wall; 




Fig. 5. Contours of streamlines in a meridional plane. Left: case (d) : (A = 1/3, Re = 2 000); Right: 
case (e) : (A = 3,Re = 2 000). The 30 contours are non-uniformly spaced for visualization purposes, 
20 equally-spaced negative contours and 10 equally-spaced positive contours for (d) and (e). 

The very different flow patterns developed by both cases (d) and (e) are highlighted by the 
very distinctive streamlines shown in Fig. O The shallow case (d) yields two vortex breakdown 
bubbles, which are off the cylinder axis while remaining attached to the free surface. The recircu- 
lation is more intense in the largest bubble, which is elongated enough to produce a recirculation 
of the fluid from the free surface all the way down to the rotating disk, and so forth. Conse- 
quently, this elongated bubble completely separate the outer region of the flow {r/R > 0.5) 
from the inner core, where the second less intense, vortex breakdown occurs. Comparatively, the 
recirculation in case (e) is fairly limited. A vortex breakdown still occurs in this case, leading to 
the formation of a small on-axis bubble, which is detached from the free surface. One can notice 
from the bending of the streamlines near the axis at the height z/H = 0.8, that another vortex 
breakdowns is in preparation — compare this effect to the similar effect on the streamlines prior 
the vortex breakdown in the closed swirling flow at Re = 1 900, see Bouffanais [69]. 

The shallow system {d) = (A = 1/3, Re = 2 000) possesses some very distinctive features as can 
be seen in Fig. [T5] (top row). The vortex lines for r/R < 0.4 being aligned with the rotation 
axis, one can easily conclude that the flow is essentially in solid-body rotation in this inner 
core region of the cavity. The meridional flow in this inner part of the cavity has a very weak 
intensity as attested by the values close to zero of the contours of the axial and radial velocity 
components — solid contour lines are positive and negative contour lines being dashed. In contrast, 
for r > 0.4 the primary recirculation of the flow is intense and predominates. The vortex lines 
bending is limited to this region and again is at the origin of the vortex breakdown appearing 
near r/R = 0.4. The boundary layer on the rotating disk is limited to the region r/R > 0.5 
and the internal jet-like shear layer close the tubular side- wall has a structure quite different 
from the cases with A = 1. Indeed in this shear layer, the axial velocity is relatively intense all 
along the tubular side- wall, unlike cases (a), {h) and (e), where the axial velocity Uz decreases 
rapidly with z/H. This observation is easily explained by the shorter height in the case studied, 
but has several considerable consequences on the flow itself. A more intense wall-jet implies a 
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more intense angular momentum jet at the free surface, which facilitates the vortex breakdown 
phenomena. The intense angular momentum free-surface jet produces an elongated recirculation 
bubble located as seen earlier, around r / R = 0.4. In this elongated bubble, the axial velocity field 
is globally positive, thereby producing an effect similar to the jet-like shear layer near the tubular 
side-wall. In turn it generates a secondary angular momentum free-surface wall jet responsible 
for the second vortex breakdown. 

Increasing the height-to-radius aspect ratio to A = 3, modifies considerably the flow dynamics 
as can be seen in Fig{T5] (bottom row). It seems clear from the previous analysis for the shallow 
case (d), that the influence of the free surface on the flow is more important when A is small. 
The proximity between the driving disk, which generates the primary flow and the free surface 
with its specific boundary conditions, leads to the complex flow dynamics earlier explained. 
Conversely, for large values of A the important distance between the spinning disk and the 
free surface is so important that it significantly reduces the effect of the presence of the free 
surface. The flow pattern presents in fine a structure very similar to the flow pattern observed 
in the closed cylinder swirling flow, except very close to the free surface. As mentioned earlier 
the recirculation bubble itself is fairly small and located on the cylinder axis likewise in the 
closed cylinder case. Finally, it is worth adding that the region of solid-body rotation is almost 
completely eliminated. Even the closest-to-the-axis vortex lines present some bending. 

As a brief conclusion of the previous study of the steady and laminar free surface swirling flows, 
it appears that the choice A = 1 for the height-to-radius aspect ratio of the cavity in presence 
of a free surface, ensures us to deal with a complex flow dynamics. Different mechanisms are 
in competition in different regions of the cavity, and in the end make the cases with A = 1 
physically more challenging and more interesting. This conclusion — valid in presence of a free 
surface — stops being valid for the closed cylinder swirling flow, and thereby explains the focus 
in the literature on cases with A > 2. 

4-3.2 Unsteady flow 

For sufficiently small Reynolds number and irrespective of A, the basic flow state is stable. As 
noted by Lopez et al. [60], when Re is increased, the basic flow state loses stability via a variety 
of Hopf bifurcations. It is worth noting that when Re tends to infinity, the stream surfaces 
and vortex surfaces — giving the streamlines and vortex lines by intersection with a meridional 
plane — must coincide. At this point, the presence of a flat free surface poses problem because of 
the constraint of having orthogonal streamlines and vortex lines on it. This apparent paradox is 
unraveled by simply letting the free surface move. Nevertheless, we know from the experiments 
carried out by Spohn et al. [8,22], that even at a Re = 6 000 the tangential flow is extremely 
intense compared to the normal one, leading to small free-surface deformations. It is very likely 
that these small amplitude deformations are not sufficient to solve our apparent paradox. At low 
Reynolds number, like those of cases (a), (6), (rf), and (e), the viscosity acts on the velocity field 
to allow the latter condition of orthogonality to be fulfilled. But when the Reynolds number is 
increased, the action of viscosity and the limited deformation of the free surface are not sufficient 
to bring back the orthogonality of the two sets of lines. Therefore, the flow must either lose its 
axisymmetry or become unsteady in order to allow to drop the orthogonality condition. The 
experiments by Spohn et al. [8] suggest that the open swirling flow first go through the unsteady 
path. 



18 




Fig. 6. Contours of streamlines in two orthogonal meridional planes, case A = 1 and Re = 6 000. Top 
row: instantaneous flow; Bottom row: mean flow. Left column: meridional plane y/R = 0; Right column: 
meridional plane x/R = 0. The 30 contours are non-uniformly spaced for visualization purposes, 20 
equally-spaced negative contours and 10 equally-spaced positive contours. 

In this section, the study is focused on the unsteady swirling flow corresponding to case (c) = 
(A = l,Re = 6 000). To our knowledge, such transitional regime at this relatively high Re 
has never been investigated nor reported in the literature. At this Reynolds number the loss 
of axisymmetry in this flat-free-surface case is evident from the observation of the contours of 
streamlines in Fig. [6] (top row). As the flow is unsteady, these recirculation bubbles are instanta- 
neous and correspond to a flow sample taken in the statistically-steady regime for t > 600 in ^ 
units, see Fig. El The streamlines of this flow sample are represented in two orthogonal meridian 
planes corresponding to y/R = and x/R = 0, in Fig. [6] (top row). Once again, the loss of 
axisymmetry appears clearly from the complex and nonaxisymmetric structure of the recircula- 
tion bubbles. Compared to the laminar and steady cases (a) and (b), the recirculation bubbles 
have their own dynamics and evolution. In a common approach to such unsteady problems, this 
complex dynamics is analyzed by the means of an averaging process, which is supplemented with 
an analysis of instantaneous flow samples equally-spaced in time. The mean flow is obtained by 
averaging 500 flow samples corresponding to successive flow states extracted every 0.25 times 
units (or equivalently every 100 iterations). Subsequently, the root-mean-square (rms) fluctua- 
tions of flow fields are calculated using the same extracted flow samples and the mean flow field 
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obtained earlier. 



The streamlines associated with the mean flow are shown in Fig. [6] (bottom row). The streamlines 
of the mean flow reveal the existence of a toroidal recirculation bubble, located off the cylinder 
axis and more surprisingly detached from the free surface. The toroidal shape and off-axis location 
of the mean recirculation bubble is in agreement with the increased-Re trend observed with cases 
(a) and (b) in Sec. 14.3. 1[ Regarding the detachment from the free surface of the mean bubble, 
it is more relevant here to notice that the instantaneous bubbles are still attached to the free 
surface. More precisely, one may notice two points: 

- in the meridional plane y/R = 0, a. small recirculation zone appears attached to the rotating 
disk for r/i? ~ 0.2; 

- in the meridional plane x/R = 0, the recirculation bubble is stretched from the free surface 
z/H = 1 down to z/H = 0.15, in a radial position r/R c::^ 0.4. 

These two observations remind the streamline patterns described in the case (d) = (A = 
1/3, Re = 2 000), with a long bubble stretching from the free surface down to the driving disk. 
The previous analysis is further conflrmed by the contours of the three velocity components 
and of the axial angular momentum for both an instantaneous flow sample and the mean flow, 
presented in Fig. [161 in the meridional plane x/ R = 0. A careful analysis of the vortex lines for 
the instantaneous flow sample shows a bending in the whole meridian plane. This bending is 
very signiflcant in the region 0.3 < r/R < 0.8 and 0.4 < z/H < 1, which corresponds to the 
limit between the primary recirculation of the flow and the secondary recirculation bubble. On 
the contrary, the vortex lines structure of the mean flow is as expected much more regular. The 
inner core region of the flow r/R < 0.4 displays a state of solid-body rotation. For both the 
instantaneous and mean flow, the jet-like shear layer along the tubular side-wall is turned into 
the interior of the flow by the free surface. Compared to the previous cases (a), (6), (d), (e), and 
also the closed swirling flow (A = l,Re = 6 000) in Bouffanais [69], the structure of this shear 
layer at Re = 6 000 reveals the presence of an intense radial jet of angular momentum at the 
free surface. 



The fluctuations of the flow with respect to its mean state have been calculated with the same 
flow samples as before. It should be noted that the fluctuation level corresponds to less than 5% 
of the maximal intensity of the respective mean flow flelds. Despite the relatively low level of 
fluctuation encountered, these fluctuations are very localized in space as can be seen in Fig. [T71 
Similarly to the mean flow flelds, the rms fluctuations of the velocity fleld and of F appear to 
be slightly nonaxisymmetric. All the three velocity components present a noticeable level of 
fluctuation near the free surface for radii close to 0.4. In this region, the free-surface radial jet of 
angular momentum reaches the inner flow, which is solid-body rotation. These fluctuations are 
therefore located in the vicinity of the stagnation point where the vortex breakdown is initiated. 
The rms-fluct nations of the vortex lines, i.e. F, are the highest in the corner region between the 
free surface and the tubular side-wall. It is in this corner, where the shear layer is turned into 
the interior by the presence of the free surface. 
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4-4 One- dimensional momentum budgets 



This section is devoted to the careful analysis of the momentum balance for the radial, azimuthal 
and axial components. This study is performed along different radial and axial lines within the 
cavity. It is of interest to determine the predominant physical terms, which are responsible 
for the complex flow dynamics depicted in the previous sections. As noticed in these previous 
sections, the structure of the flow in the inner core region is far different from the one close to 
the tubular side-wall. Similarly, the flow above the rotating driving disk has properties, which 
are not comparable to the ones below the fiat free surface. For the sake of conciseness, this 
momentum balance analysis is limited to cases (a) — (c) for which A = 1. 



4-4-1 General considerations 



The numerical integration of the Navier-Stokes equations using the spectral element method as 
described in Sec. |3l is performed in Cartesian coordinates (x, y, z) for the velocity components 
(m, f , w). Nevertheless, the axisymmetric nature of the container and of the boundary conditions 
imposed to the flow suggests the use of cylindrical coordinates. Indeed, the different physical 
terms involved in the momentum equation represented here by the Navier-Stokes equations — 
nonlinear advective term, viscous strain, pressure gradient, etc. — are better apprehended when 
expressed in cylindrical coordinates. Accordingly, all vectors and physical terms are recast as 
functions of (r, 6*,^), and for instance the velocity components are {ur,ug,Uz). 

The complete expression of the momentum equations in cylindrical coordinates reads 
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where successively appears, the velocity time derivative, the nonlinear advective term, the pres- 
sure gradient and the viscous strain. The central objective of this study is to compare the relative 
importance of some of these terms along different lines. Equation (1221) (resp. ( l23i) ) represents 
the momentum balance in the radial (resp. azimuthal) direction, and is analyzed along four 
radial lines at four different heights z/H = 0.03, 0.64, 0.95, 1, ranging from right above the 
rotating disk up to the free surface. Equation ^2M represents the momentum balance in the 
axial direction, and is analyzed along three different axial vertical lines at three radial positions 
r/R = 0.08, 0.48, 0.98, ranging from near the cylinder axis to near the tubular side-wall. 

For the sake of simplicity, some of the terms appearing in fl22l) - fl24l) are identified and denoted 
specifically in Table [51 In the sequel, the various graphs reporting the variations of these terms 
will use this nomenclature. 
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Name 



Expression 




Table 5 

Name definitions of different terms appearing in tlie momentum budget equations ()22I )-(|2 
4.4 -S Steady flows 



As discussed in Sec. I4.3.H the steady swirling flows are also fully axisymmetric. The first-order 
time derivative of the velocity fields is identically zero in Eqs. f l22l) - fl2^ . The axisymmetric 
property makes the velocity field independent of the azimuthal angle 6, and consequently all 
partial derivatives with respect to this variable vanish. Therefore, the balance in the momentum 
equation solely involves the terms described in Table [51 

As a first step, the momentum balance in the radial direction is presented in Fig. [7] for case 
(a) (left column) and case (6) (right column). At the free surface z/H = 1 (top row), the 
viscous terms are insignificant, and the flow is driven by the radial pressure deceleration, which 
is mainly counterbalanced by the centrifugal acceleration NLr2 = Ug/r and to a certain extent 
by NLrl. This analysis at the free surface still holds below the free surface at z/H = 0.94 and 
at z/H = 0.64. At this latter height, both the radial pressure deceleration and the centrifugal 
acceleration NLr2 have a lower magnitude than at the free surface but their magnitude is less 
localized than at the free surface. As expected, above the disk, at z/H = 0.03, all the terms 
have a higher magnitude and the momentum balance is more complex as only one single term 
NLr3 does not really contribute to the balance. The centrifugal acceleration NLr2 keeps its 
predominant position, but its maximum is now shifted towards the outer radial region, which 
corresponds to the region of highest angular momentum 0.8 < r/R < 1. The other acceleration 
term NLrl becomes relatively important. Both of these acceleration terms are counterbalanced 
by the radial pressure gradient and now also by the viscous strain. The importance of the viscous 
strain at this height z/H = 0.03 can easily be understood, as we are located in the viscous layer 
generated by the motion of the disk. One may add that in the inner core region of the flow, 
say r/R < 0.15, the linear trend observed for both the centrifugal acceleration and the pressure 
deceleration are well-know features of a flow in solid-body rotation, as observed with the vertical 
vortex lines in Fig. [H] (extreme right column). 



The evolution of those momentum balances along radial lines at different heights does not really 
change when increasing Re from 900 up to 1 500. But some noticeable trends are observable. For 
instance, even if the viscous strain does not play a central role at the free surface, it is worth 
noting that its effect is increased with Re while conversely it is decreased when getting closer 
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to the rotating bottom end-wall. 



As a second step, we aim at analyzing the momentum balance in the azimuthal direction along 
radial lines at the same different heights as before. The results are presented in Fig. [8] for case 
(a) (left column) and case (6) (right column). It is important to note at this point that the two 
terms NLtl and VStl involve partial derivatives with respect to the radial variable r. Given the 
fact that our solution is continuous and first-order differentiable within a spectral element and 
only continuous at the element edges, one expects some slight unphysical deformations of the 
plots associated with these two terms. Along the radial lines of interest, the spectral element 
edges are located at r/R = 0.2, 0.4, 0.6, 0.8, 0.97. As a consequence, some rapid variations of 
the terms VStl and NLtl are going to be simply disregarded in the coming discussions. 

A rapid glance at all the plots in Fig. [8] allows to conclude that the importance of the azimuthal 
momentum transfers resides in the near bottom end-wall region. The magnitude of all terms 
is over ten times smaller at z/H = 0.64, 0.95, 1, compared to z/H = 0.03. Given the solid- 
body rotation in the inner core region of the flow r/R < 0.15, most of the terms are vanishing 
small — excluding the unphysical values of VStl. 

At z/H = 0.03, one can notice the vigorous action of the viscous strain term VSt2 which literally 
drives the fluid in the viscous Ekman layer. This driving viscous term is being compensated by 
the convective terms NLtl and NLt3, and by the Coriolis term NLt2. As one gets closer to the 
corner between the rotating disk and the tubular side-wall, say r/R> 0.8, the interplay between 
the various terms is being reversed. The term VSt2, which is driving the fluid in the inner region 
of the cavity is now a dissipative term in the jet-like shear layer. Conversely, the convective term 
NLtl becomes large and is driving the fluid in the shear layer. Very close to the tubular side-wall, 
this term starts being counterbalanced by the second viscous term VStl. 

The viscous driving effect of the term VSt2 becomes insignificant at z/H = 0.64, but when 
z/H is increased, VSt2 starts growing again to reach a local maximum value at the free surface, 
but with a magnitude slightly smaller than the two other nonlinear convective terms NLtl and 
NLt2. Close to the free surface and at the free surface, the flow is primarily driven by the Coriolis 
term NLt2 together with the viscous term VSt2. Their global action is counterbalanced by the 
nonlinear convective term NLtl. This momentum balance at the free surface in the outer region 
r/R > 0.4 reflects the central effect of the free-surface jet of angular momentum. Again the 
viscous effects are more intense at the free surface when Re is increased from 900 to 1 500, 
despite the presence of the kinematic viscosity term 1/Re in their definitions. 

As a last step for these two steady flows (a) and (6), we investigate the momentum transfer in 
the axial direction, but now along three different vertical lines corresponding to three different 
radii r/R = 0.08, 0.48, 0.98. The graphs for the five different terms involved are reported in 
Fig. O We start from the top row, which is associated with the results for r/R = 0.98 in the 
jet-like shear layer and which present the highest magnitudes of all radii considered. The lack 
of axial momentum transfers is clearly visible in the region z/H > 0.4 and even inexistent 
at the free surface. Conversely, in the corner between the rotating bottom end-wall and the 
tubular side- wall the flow is driven by the axial pressure gradient, and to some extent by the two 
nonlinear terms NLal and NLa2 independently. The viscous term VSal is primarily responsible 
for counterbalancing the driving pressure effects. The second viscous term VSa2 has a less 
important effect in terms of magnitude but is driving very close to the disk (viscous Ekman 
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Fig. 7. Momentum balance in radial direction plotted along horizontal radial lines at four different 
vertical positions. First row: z/H = 1; second row: z/H = 0.95; third row: z/H = 0.64; fourth row: 
z/H = 0.03. Left column: case Re = 900; Right column: case Re = 1 500. Case A = 1. The terminology 
refers to Tab. [H 



layer), when z/H is increased it becomes dissipative as we are out of the Ekman layer but 
within the wall-jet shear layer. 

Near the medium radial position r/R= 0.48, the axial momentum transfers near the free surface 
show the relative importance of the two nonlinear terms NLal and NLa2, which are driving the 
fluid against the axial pressure gradient, which is negative as a consequence of the kinematic 
boundary condition imposing a vanishing axial velocity component at the free surface. 

Near the cylinder axis, for r/R= 0.08, the axial momentum transfers are limited and relatively 
simples. Nonlinear convective effects seem almost insignificant and the axial component of the 
flow is obtained from the balance between the axial pressure gradient and the two viscous terms: 
VSa2 in the bottom of the cavity and VSal in the top of the cavity. 
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Fig. 8. Momentum balance in azimuthal direction plotted along horizontal radial lines at four different 
vertical positions. First row: z/H = 1; second row: z/H = 0.95; third row: z/H = 0.64; fourth row: 
z/H = 0.03. Left column: case Re = 900; Right column: case Re = 1 500. Case A = 1. The terminology 
refers to Tab. [5j 



^.^.5 Unsteady flow 

The momentum balance analysis developed in the previous section for the two steady flows 
(a) and (6) cannot be directly transposed to the unsteady case (c). The two reasons for that 
are first the unsteady character requiring to account for the first-order time derivative du/dt. 
The second reason is the loss of axisymmetry of case (c) imposing to account for all the terms 
involving a partial derivative with respect to 6 in Eqs. (l22l)-(!24l). One can overcome the issue 
associated with du/dt by performing the analysis on the mean flow, which is obtained from 
the statistically steady regime. On the other hand, the momentum equations for the mean flow 
involve the Reynolds stress terms, expressing the influence of the fluctuating velocity field onto 
the dynamics of the mean field. Despite all these considerations and issues, we have deliberately 
omitted the terms involving derivatives with respect to the time and to the azimuthal coordinate, 
and we have calculated the values of the terms in Table on the same lines as in Sec. 14.4.21 By 
doing so, the objective is not to reproduce a similar analysis as with the steady cases, but more 
to investigate the evolution of the different terms for the instantaneous flow and the mean flow 
as compared to the laminar cases (a) and (6). 
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Fig. 9. Momentum balance in axial direction plotted along vertical lines at three different radial po- 
sitions. First row: r/R = 0.98; second row: r/R = 0.48; third row: r/R = 0.08. Left column: case 
Re = 900; Right column: case Re = 1 500. Case A = 1. The terminology refers to Tab.[5j 

We compare the radial terms for the instantaneous and mean flows presented in Fig. [TOl to 
their laminar and steady counterparts in Fig. [71 In terms of magnitude, the leading terms have a 
slightly higher magnitude at Re = 6 000. The general observations given in Sec. l4.4.2] remain valid 
here for z/ H = 0.64, 0.95, 1. However, the influence of the intense recirculation bubble modifies 
locally and significantly the terms in the region 0.2 < r/R < 0.4, for the instantaneous flow. 
Indeed, in this interval around the radial position r/R = 0.3 and close to the free surface, the 
centrifugal acceleration NLr2 presents a local minimum, while the decelerating radial pressure 
gradient presents a local maximum. These localized effects are directly related to the presence 
of the recirculation bubble as can be seen in Fig. [61 A similar observation can be done for the 
mean flow but the effect is much less visible. 



Close to the rotating disk, at z/H = 0.03, the variations of the five axial terms are notably 
different from their laminar counterparts, but extremely similar for the instantaneous and mean 
flows. More precisely, the viscous term VSr has mainly a dissipative action in the laminar regime, 
while it is slightly driving the flow at Re = 6 000, except very close to the tubular side- wall where 
it gets back its dissipative action in the jet-like shear layer. Moreover, the two leading terms, 
namely the centrifugal acceleration NLr2 and the radial pressure gradient, both presents a local 
maximum in the region 0.2 < r/R < 0.3, for the case (c), while it keeps growing in the laminar 
regime. This particular observation is again related to the presence of the recirculation bubbles 
in this region, which locally strongly modifles the momentum transfers. Finally, the third radial 
convective term NLr3, which is very small in the laminar regime, acquires a magnitude as 
important as the two other convective terms for r/R> 0.9. 

Let us consider now the radial variations of the flve azimuthal terms as shown in Fig. [TTl A rapid 
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overlook of all variations for the instantaneous flow sample (left column) allows to conclude to a 
general agreement with the results obtained in the laminar cases (a) and (6). The variations of 
the different terms are similar for the instantaneous and mean flows in the outer radial region, 
which implies again a relative steadiness of those terms for r/R > 0.7. On the other hand, 
the mean flow yields vanishingly small terms in the inner core region r/R < 0.4, where the 
instantaneous flow have the three nonlinear convective terms NLtl-NLt3 with a relatively high 
magnitude. The unsteady activity of those three nonlinear terms and the intense fluctuating 
activity generated by them is further discussed in Sec. 14.51 

As a last step, we compare the axial terms for the instantaneous and mean flows presented in 
Fig. [12] to their laminar and steady counterparts in Fig. [91 We start from the outer radial line 
r/R = 0.98, where the results for the instantaneous flow and the mean flow are extremely close, 
revealing an almost steady behavior of the jet-like wall shear layer surrounding the tubular side- 
wall. The comparison of these results with those of cases (a) and (b) leads to several comments. 
First, the variations of all the terms are limited to a smaller zone above the disk at Re = 
6 000. The magnitude of the axial pressure gradient is increased with Re. More surprisingly the 
axial pressure gradient is no longer counterbalanced by the convective term NLal, but is now 
counterbalanced by NLa2. By extension, one can infer that the jet-like shear layer is dominated 
by axial effects at high Reynolds number. 

For the two other radii r/R = 0.08, 0.48, the instantaneous flow terms are far different from 
their mean counterparts. As consequence, the flow in the inner core region appears much more 
unsteady than the flow in the outer region of the cavity. In addition, given the high level of 
fluctuation in the inner core region of the flow — see Fig. [17] — it appears irrelevant to further 
analyze the results for the mean flow. On the other hand, variations of the different axial terms 
for the instantaneous flow reveals that the axial momentum transfers are more important at 
high Re and are predominant in the top half of the cavity, including below the free surface. 
Also, not shown here, the viscous terms are still insigniflcant and the two axial convective terms 
NLal and NLa2, and the axial pressure gradient dominate the transfers with other unsteady 
and nonaxisymmetric terms. 

4-5 Nonaxisymmetric modes in the unsteady transitional flow 

All the previous results dealing with the transitional case (c) reveals a complex flow dynamics 
due to instabilities developing from a steady stable flow similar to the steady laminar cases (a) 
and (6). The objective of this section is to discuss the physical origin of this instability based on 
the results presented earlier and also to characterize, at least qualitatively, its effect on the flow 
fleld. 

As mentioned on several occasions in Sec. 14.31 and Sec. 14. 4[ the inner core region of the flow 
r/R < 0.2 is globally governed by a quasi-solid-body rotation and in the outer radial region, a 
wall-jet along the tubular side- wall drives the flow. As already discussed, this latter axial wall-jet 
is turned into a radial surface jet by the presence of the free surface. It seems therefore legitimate 
to consider the interfacial zone — denoted T in the sequel — between the inner core region and the 
radially- inward coming jet as prone to developing instabilities. Given the stress-free condition 
imposed on the free surface, the effect of the instabilities should persist all the way to the free 
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Fig. 10. Momentum balance in radial direction plotted along horizontal radial lines at four different 
vertical positions. First row: z/H = 1; second row: z/H = 0.95; third row: z/H = 0.64; fourth row: 
z/H = 0.03. Left column: instantaneous flow; Right column: mean flow. Case A = 1 and Re = 6000. 
The terminology refers to Tab. [5j 



surface itself. Consequently, we consider the variations at the free surface z/H = 1 of the radial 
and azimuthal velocity components — the axial component vanishes at the free surface because 
of the kinematic boundary condition on it — and of the axial angular momentum F = rug. These 
variations for the instantaneous flow are shown in Fig. [18] (top row), while the corresponding 
rms fluctuations are presented in the row below. 

The most significative feature highlighted in these graphs is the presence of an annular region 
with 0.3 < r/R < 0.4, where the fluctuations of Ur, ug and F are intense. Such intense fluctuating 
activity brings the interfacial zone T to light. Focusing now on the rms fluctuations of the axial 
angular momentum, one may notice the presence of a second outer annular region comprised in 
the interval 0.6 < r/R < 0.9, which is nonaxisymmetric and relatively intense. This outer region 
of intense rms fluctuations for F correspond to the zone where the jet-like shear layer is turned 
into the interior by the free surface. This outer wall jet injects high-angular-momentum fluid 
towards the cylinder axis. This radial jet impinges on the cylindrical core of the fluid that is in 
solid-body rotation. When increasing the Reynolds number, the radius of the cylindrical core in 
solid-body rotation is reduced, while the intensity of the radial jet is increased. Above a given 
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Fig. 11. Momentum balance in azimuthal direction plotted along horizontal radial lines at four different 
vertical positions. First row: zjH = 1; second row: z/H = 0.95; third row: z/H = 0.64; fourth row: 
z/H = 0.03. Left column: instantaneous flow; Right column: mean flow. Case A = 1 and Re = 6000. 
The terminology refers to Tab. [5j 



value of the Reynolds number, the action of the impingement of the radial jet on the inner core 
region starts developing unstable modes. 

The origin of these unstable modes is to be found in the analysis of the momentum transfers 
performed in the previous section. Returning on the radial variations of the five azimuthal 
terms below the free surface at z/H = 0.95, Fig. [13] reproduces these variations for Re = 
900, 1500, 6 000 (instantaneous flow), 6 000 (mean flow) from top to bottom. All these results 
have been shown separately before and are now shown together to facilitate the discussion. 
The flow in the outer region r/R > 0.5 has a strong steady character given the fact that the 
variations for the mean flow are fairly close to those of the instantaneous flow. The Coriolis term 
NLt2 = —UrUg/r keeps the same radially-outward decreasing trend. On the other hand, the 
two other nonlinear terms NLtl and NLtS develop opposed and equally-intense peaks around 
r/R = 0.9. These opposite peaks have an increasing intensity with the Reynolds number. In the 
interfacial zone T, all the terms involved present brutal variations and changes of behavior, which 
give another characterization of this interfacial zone T. In the inner core region r/R < 0.4, the 
flow possesses a strong unsteady character brought to light by the vanishingly small values of the 
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Fig. 12. Momentum balance in axial direction plotted along vertical lines at three different radial 
positions. First row: r/R = 0.98; second row: r/R = 0.48; third row: r/R = 0.08. Left column: 
instantaneous flow; Right column: mean flow. Case A = 1 and Re = 6 000. The terminology refers to 
Tab. El 

various terms for the mean flow. Consequently, the unstable azimuthal modes are to be found into 
this inner cylindrical region. Indeed, one may notice that the convective term NLtl = —Urdue/dr 
and the Coriolis term NLt2 = —UrUg/r are negative and have a low magnitude in the laminar 
cases (a) and (b), but acquires large positive values at Re = 6 000. Moreover, all the three 
other terms NLt3, VStl and VSt2 have very low magnitudes, and thus cannot counterbalance 
the azimuthal momentum injected by NLtl and the Coriolis term NLt2. Only an unsteady and 
nonaxisymmetric flow can support such azimuthal momentum effects. The effect of these two 
destabilizing terms on the flow apparently leads to the formation of azimuthal rotating waves 
superimposed to the stable base flow. The variations at the free surface z/H = 1, of the radial 
and azimuthal velocity components, and of the axial angular momentum shown in Fig. [18] (top 
row), suggest the conjugate effect of several rotating waves. These rotating waves correspond to 
even azimuthal Fourier modes, mainly n = 2 and n = 4. 



5 Conclusions 

The incompressible flow of a viscous fluid enclosed in a cylindrical container with an open top flat 
surface and driven by the constant rotation of the bottom wall has been thoroughly investigated. 
The top surface of the cylindrical cavity is left open with a stress- free boundary condition imposed 
on it. No-slip condition is imposed on the side- wall and also on the rotating bottom end- wall by 
means of a regularized angular velocity proflle. More speciflcally, the stress-free top surface is 
maintained flxed and flat. 
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Fig. 13. Momentum terms in azimuthal direction plotted along the horizontal radial line at z/H = 0.95. 
Prom top to bottom: case (a); case (6); case (c) instantaneous; case (c) mean flow. The terminology 
refers to Tab. [3 

New flow states have been investigated based on a fully three-dimensional solution of the Navier- 
Stokes equations for the free-surface cylindrical swirling flow, without resorting to any symmetry 
property unlike all other results available in the literature. To our knowledge, the present study 
delivers the most general available results for this flat-free-surface problem due to its original 
mathematical treatment. 



Five different cases corresponding to different pairs of governing parameters (A, Re) have been 
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considered. The Reynolds regime corresponds to transitional flows with some incursions in the 
fully laminar regime. Both steady and unsteady non-oscillatory swirling flows are considered 
with a particular emphasis on the case (A = 1, Re = 6 000). Of great concern to this study is the 
question of space resolution. This is particularly important for the bifurcated case at Re = 6 000. 
Convergence tests in space and time have been carried out on this upmost problematic case, and 
optimal values of the polynomial degree and time-step have been deduced. 

The evolution of the total kinetic energy of this open flow has been carefully studied for increasing 
Reynolds numbers and has been compared to the results for the closed swirling flow. The presence 
of the free surface on the top of the cylinder is found to strongly modify the observed trend: the 
total kinetic energy is increased with Re in the open cylinder case, while the converse is observed 
in the closed cylinder case. A physical analysis of the energetic action of the surface layer below 
the free surface allows to justify the above results. A comprehensive physical description of all 
flow states has been given with particular emphasis on the vortex breakdown bubbles and on the 
structure of the vortex lines. The unsteady case at Re = 6 000 has retained more attention, given 
its unsteady transitional character. The mean flow and the corresponding rms fluctuations have 
been calculated and the results analyzed accordingly. The momentum transfers in the radial, 
azimuthal and axial directions have been studied along various one-dimensional lines. For the 
transitional case at Re = 6 000, the flow in an inner cylindrical core is in solid-body rotation, 
while the outer radial layer is dominated by the jet-like shear layer along the tubular side- wall. 
This axial wall-jet is turned into a radial jet of angular momentum, which prevails all the way 
up to the free surface. The impingement of this radial jet onto the inner cylindrical core in 
solid-body rotation leads to the development of unstable azimuthal modes. The nonlinear terms, 
which includes a Coriolis effect, responsible for the development of these unstable modes have 
been found using the azimuthal momentum imbalance below the free surface. These unstable 
modes take the form of even-order azimuthal rotating waves. 
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Fig. 14. Contours in a meridional plane for the case A = 1. Top row: case Re = 900; Bottom row: case Re = 1 500. From left column to right 
column: radial velocity component Ur] axial velocity component w = Uz', azimuthal velocity component uq; axial angular momentum component 
r = rug. The 35 contours are uniformly spaced, between —0.1 and 0.145 for Ur] and between —0.08 and 0.115 for Uz', The 50 contours are 
uniformly spaced, between and 1 for ug and F. 
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Fig. 15. Contours in a meridional plane for the case (d) : (A = 1/3, Re = 2 000) (top row) and case (e) : (A = 3, Re = 2 000) (bottom row). 
From left column to right column: radial velocity component Ur', axial velocity component w = Uz] azimuthal velocity component uq] axial 
angular momentum component F = rug. All contours are uniformly spaced; 35 contours between —0.17 and 0.15 for case (d); 90 contours 
between —0.03 and 0.15 for Ur case (e); 35 contours between —0.1 and 0.12 for Uz case (d); 35 contours between —0.05 and 0.12 for Uz case (e); 
50 contours between and 1 for uq and F, cases (d) and (e). 





Fig. 16. Contours in a meridional plane for the case A = 1 and Re = 6 000. Top row: instantaneous flow; Bottom row: mean flow. From left 
column to right column: radial velocity component Ur] axial velocity component Uz] azimuthal velocity component uq; axial angular momentum 
r = ruQ. The 50 contours are uniformly spaced, between —0.13 and 0.16 for n^; between —0.09 and 0.14 for Uz\ and between and 1 for ug. 
The 100 contours are uniformly spaced between and 1 for F. 




Fig. 17. Contours in a meridional plane for the case A = 1 and Re = 6 000. Top row: in the meridian plane y/R = 0; Bottom row: in the 
meridian plane x/R = 0. From left column to right column: rms fluctuations of radial velocity component Ur] rms fluctuations of the axial 
velocity component Uz] rms fluctuations of the azimuthal velocity component ug; and rms fluctuations of the axial angular momentum F = rug. 
The 20 contours are uniformly spaced, between and 0.002 for rms-Ur! between and 0.0024 for ims-Uz] between and 0.004 for rms-u^; and 
between and 0.001 for rms-F. 



Fig. 18. Contours on the free surface z/H = 1 for the case A = 1 and Re = 6 000. Top row: instantaneous flow; Bottom row: rms fluctuations. 
Left column: radial velocity component u^; Central column: azimuthal velocity component uq] Right column: axial component of the angular 
momentum F = rug. The 15 contours are uniformly spaced, between —0.13 and 0.03 for Ur] between and 0.35 for uq] and between and 0.2 
for F. The 20 contours are uniformly spaced, between and 0.002 for rms-Ur-; between and 0.024 for rms-ug; and between and 0.001 for 
rms-F. 



